library(ggplot2)
library(ggpubr)
theme_customize=theme_bw()+
  theme(panel.border=element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(colour = "black"),
        legend.title = element_blank())+
  theme(axis.text.x=element_text(size=14,angle=0,color="Black"),
        axis.text.y=element_text(size=14,angle=0,color="Black"))+
  theme(axis.title.x = element_text(size = 15, color = "Black"),
        axis.title.y = element_text(size = 15,  color = "Black"))

multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)
  
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  
  numPlots = length(plots)
  
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  
  if (numPlots==1) {
    print(plots[[1]])
    
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

library(vegan)
setwd("D:/入侵/数据/分析2021/2 抽吸1")
dat_sim_site=read.csv("sim-each site.csv")
dat_sim_site=t(dat_sim_site)
shannon_index <- diversity(dat_sim_site, index = 'shannon', base = exp(1))
write.csv(shannon_index,file = "shannon_index.csv",row.names = T)


library("reshape2")
results=read.csv("results.csv",header=TRUE)
dat1<-melt(results,
           id.vars = c('Community','Site'),
           variable.name='ind',
           value.name='val')

library(dplyr)
list_name1 <-  unique(dat1$ind)
##normality 
sink("normal-test result.txt") 
for(i in list_name1){
  print(paste(i))
  tmp<-shapiro.test((dat1 %>% filter(ind == i))$val)
  print(tmp)
}
sink()
##homogeneity of variance 
library(car)
sink("leveneTest result.txt")  
for(i in list_name1){
  print(paste(i))
  tmp<-dat1 %>% filter(ind == i) %>% leveneTest(val~Community,data=.)
  print(tmp)
}
sink()

library(agricolae)
sink("aov result.txt") 
for(i in list_name1){
  print(paste(i))
  mod<-dat1 %>% filter(ind == i) %>% aov(val~Community,data=.)
  print(summary(mod))
  out<- LSD.test(mod,"Community",p.adj="none" )
  print(out$groups)
}
sink()

sink("Kruskal result.txt")  
for(i in list_name1){
  print(paste(i))
  mod<-dat1 %>% filter(ind == i) %>% with(.,kruskal(val,Community,group=TRUE, main=paste(i,sep = " ")))
  print(mod)
}
sink()

####Drawing figures
###Richness index
rich=read.csv("rich for bar.csv",header=TRUE)
rich$Community = factor(rich$Community, levels=c('OP','IP','IPS','IS','RP'))

p_Richness=ggplot(rich,aes(x=Community,y=Richness,fill=Community))+
  theme_customize+
  theme(legend.position="none")+
  geom_bar(position="dodge",stat="identity",width=0.5)+
  geom_text(aes(y=Richness+1.96*SE+0.6,label=groups),position=position_dodge(width = 0.9))+
  xlab("Plant Community")+
  ylab("Richness")+
  geom_errorbar(aes(ymax = Richness+1.96*SE, ymin = Richness-1.96*SE), 
                position = position_dodge(0.9), width = 0.15)


###Shannon-Wiener index
shan=read.csv("shan for bar.csv",header=TRUE)
shan$Community = factor(shan$Community, levels=c('OP','IP','IPS','IS','RP'))

p_Shannon=ggplot(shan,aes(x=Community,y=Shannon))+
  theme_customize+
  theme(legend.position = c(0.86,0.88))+
  geom_bar(position="dodge",stat="identity",width=0.5)+
  geom_text(aes(y=Shannon+1.96*SE+0.05,label=groups),position=position_dodge(width = 0.9))+
  xlab("Plant Community")+
  ylab("Shannon")+
  geom_errorbar(aes(ymax = Shannon+1.96*SE, ymin = Shannon-1.96*SE), 
                position = position_dodge(0.9), width = 0.15)


###Density
dens=read.csv("den for bar.csv",header=TRUE)
dens$Community = factor(dens$Community, levels=c('OP','IP','IPS','IS','RP'))

p_Density=ggplot(dens,aes(x=Community,y=Density))+
  theme_customize+
  theme(legend.position = c(0.86,0.88))+
  geom_bar(position="dodge",stat="identity",width=0.5)+
  geom_text(aes(y=Density+1.96*SE+20,label=groups),position=position_dodge(width = 0.9))+ 
  xlab("Plant Community")+
  ylab("Density (numbers per m^2)")+
  geom_errorbar(aes(ymax = Density+1.96*SE, ymin = Density-1.96*SE), 
                position = position_dodge(0.9), width = 0.15)


multiplot(p_Richness,p_Density,p_Shannon, cols=3)  #10*3.5